
; Define loc
inpath_b   = "../../../daily_to_monthly/"
inpath_a   = "../../../../2018_2018/daily_to_monthly/"
outpath   = "./"

; Define time

do y=2018,2018

do m=1,12

t_name_in = sprinti("%0.4i",y) + "-" + sprinti("%0.2i",m)\
		+ "-" + sprinti("%0.2i",1) + "_" + sprinti("%0.2i",0)\
		+ "_00_00.nc"

print(t_name_in)


;****************************************************************
;wind: m/s; specific humidy Kg/Kg; ps hPa; omega Pa/s; lev hPa**

; ******************Path before expansion************************
fin = addfile(inpath_b + "wrfout_monthly_d02_" + t_name_in, "r")

; ext data
tc_eta       = wrf_user_getvar(fin,"tc",-1) + 273.15
tc_eta@units = "degK"
p_eta        = wrf_user_getvar(fin,"p",-1)
wa_eta       = wrf_user_getvar(fin,"wa",-1)
w_eta   	 = w_to_omega(wa_eta, p_eta, tc_eta)
delete([/wa_eta, tc_eta, p_eta/])

q_eta   = wrf_user_getvar(fin,"QVAPOR",-1)     ; The variable to interpolate
q_eta 	 = q_eta/(q_eta+1)
q_eta@units =  "specific humidity"
uvmet_eta    = wrf_user_getvar(fin,"uvmet",-1)
u_eta   = uvmet_eta(0,:,:,:,:)
v_eta   = uvmet_eta(1,:,:,:,:)

; inter level
vert_coord          = "pres"
plev       = ispan(1000, 100, 10)		; the inter data is consistent with this order
plev@units =  "Pressure"
opts             = True
opts@extrapolate = False 		; under ground
opts@field_type  = "pres"
opts@logP        = True   ; Set true when inter press

; inter
qb      	= wrf_user_vert_interp(fin,q_eta,vert_coord,plev,opts)
wb      	= wrf_user_vert_interp(fin,w_eta,vert_coord,plev,opts)
ub   		= wrf_user_vert_interp(fin,u_eta,vert_coord,plev,opts)
vb  		= wrf_user_vert_interp(fin,v_eta,vert_coord,plev,opts)
psb    		= fin->PSFC
psb			= psb/100 ;hpa
psb@units	= "hpa"
delete([/w_eta, q_eta, u_eta, v_eta, uvmet_eta/])

; ******************Path after expansion************************
fin = addfile(inpath_a + "wrfout_monthly_d02_" + t_name_in, "r")

; ext data
tc_eta       = wrf_user_getvar(fin,"tc",-1) + 273.15
tc_eta@units = "degK"
p_eta        = wrf_user_getvar(fin,"p",-1)
wa_eta       = wrf_user_getvar(fin,"wa",-1)
w_eta   	 = w_to_omega(wa_eta, p_eta, tc_eta)
delete([/wa_eta, tc_eta, p_eta/])

q_eta   = wrf_user_getvar(fin,"QVAPOR",-1)     ; The variable to interpolate
q_eta 	 = q_eta/(q_eta+1)
q_eta@units =  "specific humidity"
uvmet_eta    = wrf_user_getvar(fin,"uvmet",-1)
u_eta   = uvmet_eta(0,:,:,:,:)
v_eta   = uvmet_eta(1,:,:,:,:)

; inter
qa      	= wrf_user_vert_interp(fin,q_eta,vert_coord,plev,opts)
wa      	= wrf_user_vert_interp(fin,w_eta,vert_coord,plev,opts)
ua   		= wrf_user_vert_interp(fin,u_eta,vert_coord,plev,opts)
va  		= wrf_user_vert_interp(fin,v_eta,vert_coord,plev,opts)
psa    		= fin->PSFC
psa			= psa/100 ;hpa
psa@units	= "hpa"
delete([/w_eta, q_eta, u_eta, v_eta, uvmet_eta/])

; ******************diff between after and before expansion************************
qf = qa - qb
wf = wa - wb
uf = ua - ub
vf = va - vb

copy_VarCoords(qb,qf)
copy_VarCoords(wb,wf)
copy_VarCoords(ub,uf)
copy_VarCoords(vb,vf)

; plev dx dy
plev   = qb&interp_levels
dx     = 5000
dx@units = "meter"
dy     = 5000
dy@units = "meter"
var_size = dimsizes(qb)
nx = var_size(2)
ny = var_size(3)
linlog = 1
g = 9.8
pbot = max(plev)
ptop = min(plev)

;*******************Water vapor transport flux volumn***********
; **********************before**********************
uqb = ub
vqb = vb
uqb = ub*qb
vqb = vb*qb

UQb_vint = 1./g * 100 * vibeta (plev,uqb(Time|:,south_north|:,west_east|:,interp_levels|:),linlog,psb,pbot,ptop) ; kg/(m*s) 
VQb_vint = 1./g * 100 * vibeta (plev,vqb(Time|:,south_north|:,west_east|:,interp_levels|:),linlog,psb,pbot,ptop)
copy_VarCoords(psb,UQb_vint)
copy_VarCoords(psb,VQb_vint)
UQb_vint@units = "kg m-1 s-1"
VQb_vint@units = "kg m-1 s-1"

; **********************after**********************
uqa = ua
vqa = va
uqa = ua*qa
vqa = va*qa

UQa_vint = 1./g * 100 * vibeta (plev,uqa(Time|:,south_north|:,west_east|:,interp_levels|:),linlog,psa,pbot,ptop) ; kg/(m*s) 
VQa_vint = 1./g * 100 * vibeta (plev,vqa(Time|:,south_north|:,west_east|:,interp_levels|:),linlog,psa,pbot,ptop)
copy_VarCoords(psa,UQa_vint)
copy_VarCoords(psa,VQa_vint)
UQa_vint@units = "kg m-1 s-1"
VQa_vint@units = "kg m-1 s-1"

;*******************moisture advection, U V W****************
; moisture advection before and after expansion

; **********************before**********************
; -wbdqbdp
wbdqbdp = wb;
dqbdp  = center_finite_diff_n(qb, plev, False, 0, 1)
wbdqbdp = wb * dqbdp

; -ubdqbdx 
ubdqbdx = ub ;(Time|:,south_north|:,west_east|:,interp_levels|:)
do  i = 0, nx-1
    ubdqbdx(:,:,i,:) = center_finite_diff_n(qb(:,:,i,:), dx, False, 0, 2)
end do
ubdqbdx = ub * ubdqbdx

; -vbdqbdy
vbdqbdy = vb ;(Time|:,south_north|:,west_east|:,interp_levels|:)
do i = 0,ny-1
   vbdqbdy(:,:,:,i) = center_finite_diff_n(qb(:,:,:,i), dy, False, 0, 2)
end do
vbdqbdy = vb * vbdqbdy

wbdqbdp_vint = -1./g * 86400 * vibeta (plev,wbdqbdp(Time|:,south_north|:,west_east|:,interp_levels|:),linlog,psb,pbot,ptop) ;to mm/day
ubdqbdx_vint = -1./g * 100 * 86400 * vibeta (plev,ubdqbdx(Time|:,south_north|:,west_east|:,interp_levels|:),linlog,psb,pbot,ptop) ;to mm/day
vbdqbdy_vint = -1./g * 100 * 86400 * vibeta (plev,vbdqbdy(Time|:,south_north|:,west_east|:,interp_levels|:),linlog,psb,pbot,ptop) ;to mm/day

copy_VarCoords(psb,wbdqbdp_vint)
copy_VarCoords(psb,ubdqbdx_vint)
copy_VarCoords(psb,vbdqbdy_vint)

; **********************after**********************
; -wadqadp
wadqadp = wa;
dqadp  = center_finite_diff_n(qa, plev, False, 0, 1)
wadqadp = wa * dqadp

; -uadqadx 
uadqadx = ua ;(Time|:,south_north|:,west_east|:,interp_levels|:)
do  i = 0, nx-1
    uadqadx(:,:,i,:) = center_finite_diff_n(qa(:,:,i,:), dx, False, 0, 2)
end do
uadqadx = ua * uadqadx

; -vadqady
vadqady = va ;(Time|:,south_north|:,west_east|:,interp_levels|:)
do i = 0,ny-1
   vadqady(:,:,:,i) = center_finite_diff_n(qa(:,:,:,i), dy, False, 0, 2)
end do
vadqady = va * vadqady

wadqadp_vint = -1./g * 86400 * vibeta (plev,wadqadp(Time|:,south_north|:,west_east|:,interp_levels|:),linlog,psa,pbot,ptop) ;to mm/day
uadqadx_vint = -1./g * 100 * 86400 * vibeta (plev,uadqadx(Time|:,south_north|:,west_east|:,interp_levels|:),linlog,psa,pbot,ptop) ;to mm/day
vadqady_vint = -1./g * 100 * 86400 * vibeta (plev,vadqady(Time|:,south_north|:,west_east|:,interp_levels|:),linlog,psa,pbot,ptop) ;to mm/day

copy_VarCoords(psa,wadqadp_vint)
copy_VarCoords(psa,uadqadx_vint)
copy_VarCoords(psa,vadqady_vint)

;*******************Water vapor budget volumn****************
; P = -<delt(Vq)> + E + residual  ***  -<delt(Vq)> = -wdqdp-udqdx-vdqdy

;*******************-<wdqdp>'****************
; -<wdqdp>' ~= -wbdqfdp -wfdqbdp -wfdqfdp  (qf = qa - qb; wf = wa -wb)
wbdqfdp = wb
wfdqbdp = wf
wfdqfdp = wf
dqfdp   = center_finite_diff_n(qf ,plev,False,0, 1)
dqbdp   = center_finite_diff_n(qb ,plev,False,0, 1)
wbdqfdp = wb * dqfdp
wfdqbdp = wf * dqbdp
wfdqfdp = wf * dqfdp

Ther_w = -1./g * 86400 * vibeta (plev,wbdqfdp(Time|:,south_north|:,west_east|:,interp_levels|:),linlog,psb,pbot,ptop) ;to mm/day
Dyna_w = -1./g * 86400 * vibeta (plev,wfdqbdp(Time|:,south_north|:,west_east|:,interp_levels|:),linlog,psb,pbot,ptop) ;to mm/day
Res_w  = -1./g * 86400 * vibeta (plev,wfdqfdp(Time|:,south_north|:,west_east|:,interp_levels|:),linlog,psb,pbot,ptop) ;to mm/day

copy_VarCoords(psb,Ther_w)
copy_VarCoords(psb,Dyna_w)
copy_VarCoords(psb,Res_w)

;*******************-<udqdp>'****************
; -<udqdp>' ~= -ubdqfdp -ufdqbdp -ufdqfdp  (qf = qa - qb; uf = wa -wb)

ubdqfdx = ub ;(Time|:,south_north|:,west_east|:,interp_levels|:)
do  i = 0, nx-1
    ubdqfdx(:,:,i,:) = center_finite_diff_n(qf(:,:,i,:), dx, False, 0, 2)
end do
ubdqfdx = ub * ubdqfdx 

ufdqbdx = uf 
do  i = 0, nx-1
    ufdqbdx(:,:,i,:) = center_finite_diff_n(qb(:,:,i,:), dx, False, 0, 2)
end do
ufdqbdx = uf * ufdqbdx 

ufdqfdx = uf 
do  i = 0, nx-1
    ufdqfdx(:,:,i,:) = center_finite_diff_n(qf(:,:,i,:), dx, False, 0, 2)
end do
ufdqfdx = uf * ufdqfdx 

Ther_u = -1./g * 100 * 86400 * vibeta (plev,ubdqfdx(Time|:,south_north|:,west_east|:,interp_levels|:),linlog,psb,pbot,ptop)
Dyna_u = -1./g * 100 * 86400 * vibeta (plev,ufdqbdx(Time|:,south_north|:,west_east|:,interp_levels|:),linlog,psb,pbot,ptop)
Res_u  = -1./g * 100 * 86400 * vibeta (plev,ufdqfdx(Time|:,south_north|:,west_east|:,interp_levels|:),linlog,psb,pbot,ptop) ;to mm/day

copy_VarCoords(psb,Ther_u)
copy_VarCoords(psb,Dyna_u)
copy_VarCoords(psb,Res_u)

;*******************-<vdqdp>'****************
; -<vdqdp>' ~= -vbdqfdp -vfdqbdp -vfdqfdp  (qf = qa - qb; vf = wa -wb)

vbdqfdy = vb ;(Time|:,south_north|:,west_east|:,interp_levels|:)
do  i = 0, ny-1
    vbdqfdy(:,:,:,i) = center_finite_diff_n(qf(:,:,:,i), dy, False, 0, 2)
end do
vbdqfdy = vb * vbdqfdy 

vfdqbdy = vf 
do  i = 0, ny-1
    vfdqbdy(:,:,:,i) = center_finite_diff_n(qb(:,:,:,i), dy, False, 0, 2)
end do
vfdqbdy = vf * vfdqbdy 

vfdqfdy = vf 
do  i = 0, ny-1
    vfdqfdy(:,:,:,i) = center_finite_diff_n(qf(:,:,:,i), dy, False, 0, 2)
end do
vfdqfdy = vf * vfdqfdy 

Ther_v = -1./g * 100 * 86400 * vibeta (plev,vbdqfdy(Time|:,south_north|:,west_east|:,interp_levels|:),linlog,psb,pbot,ptop)
Dyna_v = -1./g * 100 * 86400 * vibeta (plev,vfdqbdy(Time|:,south_north|:,west_east|:,interp_levels|:),linlog,psb,pbot,ptop)
Res_v  = -1./g * 100 * 86400 * vibeta (plev,vfdqfdy(Time|:,south_north|:,west_east|:,interp_levels|:),linlog,psb,pbot,ptop) ;to mm/day

copy_VarCoords(psb,Ther_v)
copy_VarCoords(psb,Dyna_v)
copy_VarCoords(psb,Res_v)

;==================================output=========================================
ubdqbdx_vint@units = "mm/day"
vbdqbdy_vint@units = "mm/day"
wbdqbdp_vint@units = "mm/day"
uadqadx_vint@units = "mm/day"
vadqady_vint@units = "mm/day"
wadqadp_vint@units = "mm/day"
Ther_u@units     = "mm/day"
Dyna_u@units     = "mm/day"
Ther_v@units     = "mm/day"
Dyna_v@units     = "mm/day"
Ther_w@units     = "mm/day"
Dyna_w@units     = "mm/day"
Res_u@units      = "mm/day"
Res_v@units      = "mm/day"
Res_w@units      = "mm/day"

out_fname = outpath + "wrf_moisture_budget_monthly_d02_" + t_name_in
system("rm -rf "+out_fname)
fout = addfile( out_fname, "c")
fout->UQb_vint  	= UQb_vint
fout->VQb_vint  	= VQb_vint
fout->UQa_vint  	= UQa_vint
fout->VQa_vint  	= VQa_vint
fout->W_before  	= wbdqbdp_vint
fout->U_before 		= ubdqbdx_vint
fout->V_before 		= vbdqbdy_vint
fout->W_after  		= wadqadp_vint
fout->U_after  		= uadqadx_vint
fout->V_after  		= vadqady_vint
fout->W_Ther  		= Ther_w
fout->W_Dyna   		= Dyna_w
fout->U_Ther  		= Ther_u
fout->U_Dyna  		= Dyna_u
fout->V_Ther  	 	= Ther_v
fout->V_Dyna  		= Dyna_v
fout->U_Res   		= Res_u
fout->V_Res  		= Res_v
fout->W_Res  		= Res_w

end do

end do











